Surface Instability of Icicles 
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Quantitatively-unexplained stationary waves or ridges often encircle icicles. Such waves form 
when roughly 0.1 mm-thick layers of water flow down the icicle. These waves typically have a 
wavelength of 1cm approximately independent of external temperature, icicle thickness, and the 
volumetric rate of water flow. In this paper we show that these waves can not be obtained by naive 
Mullins-Sekerka instability, but are caused by a quite new surface instability related to the thermal 
diffusion and hydrodynamic effect of thin water flow. 

PACS numbers: 47.20.Hw, 81.30.Fb 



I. INTRODUCTION 

Interesting wave patterns often form on growing icicles 
that are covered with a thin layer of flowing water (See 
Fig.l) . For many of these patterns, the wavelength has 
a Gaussian distribution centered at approximately 8 mm; 
however, despite their common occurrence, there is no 
quantitative explanation for this wavelength distribution 
0,0. These waves are associated with the growth of the 
icicles and the flow of fluid along the icicle. Hence, there 
are several processes occurring that should be consid- 
ered. These include crystallization from the melt, latent 
heating at the ice-melt interface, laminar flow with two 
interfaces (ice- melt and melt-air), evaporation of liquid, 
and transport of heat through the surrounding air. That 
the waves tend to encircle the icicle clearly indicates the 
importance of gravity-induced flow, although the specific 
interactions between the flow, the ice growth, and the 
heat flow through both interfaces must be considered. 




FIG. 1: Waves on icicle. 
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In the study of crystal growth, such a surface instabil- 
ity is usually explained by Mullins-Sekerka theory (here- 
after abbreviated as MS). MS theory is based on two 
observations: Laplace instability and Gibbs-Thomson ef- 
fect. (For detailed explanation, refer to text books and 
original paper by Mullins and Sekerka cited in ||.) 

To a good approximation, the ice in the icicle has a 
uniform temperature of 273 K, thus temperature gradi- 
ents into the ice are insignificant, and the external tem- 
perature field is time independent and satisfies Laplace's 
equation. We further assume that the external tempera- 
ture is below 273 K; i.e., the ice is not melting on average. 
At a convex point, the temperature gradient is higher. 
Because heat flow is proportional to the gradient of tem- 
perature, the larger heat flow at a convex point rapidly 
removes latent heat from the ice surface thus allowing 
the convex points to increase in size rapidly. Conversely, 
concave regions grow relatively slowly. This phenomenon 
suggests that short-wavelength fluctuations increase in 
amplitudes more rapidly than long wavelength fluctua- 
tions. We refer to this as the Laplace instability. 

Next we explain the Gibbs-Thomson effect. The sur- 
face of a solid object has its own energy per area called 
surface free energy. If a molecule attaches itself to the 
surface near a convex point, the surface area increases, 
resulting in an increase in energy. On the other hand, if a 
molecule becomes attached to the surface near a concave 
point, absorption of the molecule makes the surface 
area smaller. Therefore, absorption of molecule at a 
concave point is more energy-efficient than is absorption 
at a convex point. For this reason, the melting point 
depends on the curvature of an object, i.e., the shape of 
the surface area. The melting point is lower at a convex 
surface (easy to melt) and is higher at a concave surface 
(hard to melt). Such an effect suppresses the fluctuation 
and makes surface flat. This is called Gibbs-Thomson 
effect, which is opposite to that of Laplace instability. 
Laplace instability enhances shorter wavelength fluc- 
tuation, and Gibbs-Thomson effect suppresses shorter 
wavelength fluctuation. From these two effects, we 
have fluctuation of specific wavelength mainly. These 
two effects are incorporated in Mullins-Sekerka's theory 
|§, B. However, the simple application of this theory 
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cannot take place in the case of icicles. The reasons 
are as follows.: Firstly, the water layers on icicles are 
too thin to produce a Laplace instability. Because the 
instability requires that the water thickness be larger 
than the wavelength of the fluctuation. Secondly, the 
curvature is too small. (Wavelength 1 cm is much larger 
than 10 [im which GT effect requires to be effective.) 
Therefore the nature of wave patterns along icicles can 
not be explained by naive MS theory. 

We neglect the fluid instabilities that require turbu- 
lence because there should be no turbulence in these wa- 
ter films H . This is because the layers are only about 100 
fxni thick and the flow speeds are about 2-4cm/s; the re- 
sulting Reynolds numbers are only about 1 and the flow 
is laminar. The hydrodynamics in thin water layer is also 
discussed by Wettlaufer et.al to explain the premelting 
dynamics, but the discussion here is essentially different 

For our analysis, we assume water flow on a ramp (See 
Fig. 2) because it is simpler to treat and relevant exper- 
iments for this geometry are available Much of the 
same processes and relative length scales occur in both 
system because the water-layer thickness (~ 10~ 4 m) is 
much smaller than the radius of icicle (~ 10 _2 m). Fur- 
thermore, Matsuda [Q, observed wave patterns on such 
an ice ramp; for example, at 8 = n/2, the wavelength 
was about 8 mm. 




FIG. 2: Flow on a ramp inclined at 6 degrees. 

Liquid flow of a thin water on a flat ramp(Benney's liq- 
uid film) produces waves 0(see also ||); however, these 
waves travel down the ramp and thus are unlike the case 
on icicles. Due to the explicit calculation, the wavelength 
of these travelling surface waves is about 1 cm, which 
agrees with wavelength along icicles, but they move at 
about 4-8 cm/s which is twice the speed of the fluid. 
Benney's wave is caused by gravity and surface tension, 
but it is unclear how it applies to the standing waves 
on ice unless the travelling waves can become pinned to 
a fixed location; such a pinning mechanism has not yet 
been proposed. 

Our approach is to assume static flow with small rip- 
ples on the ramp surface, and then calculate the growth 
rate for the ripples by solving the thermal diffusion equa- 
tion in the background fluid. In section 2, we discuss the 
fluid dynamics of a thin water flowing along a ramp, and 




FIG. 3: Two boundaries: SL and AL. 



then in section 3, we couple the thermal diffusion pro- 
cess to the flow. The thermal diffusion in air is solved in 
section 4, and in section 5 we discuss the growth rate of 
fluctuation on icicles by combining the solutions for ther- 
mal diffusion equations in two regions; air and water. 



II. HYDRODYNAMICS IN A THIN-LAYER OF 
WATER 

We consider the fluid mechanics of a thin water layer 
with depth h(x) as sketched in Fig. 3. Over each wave- 
length, the average depth is ho- There are two material 
boundaries, one is the solid-liquid boundary SL, and the 
other is the air-liquid boundary AL. The x-axis is along 
the ramp and increases in the downhill direction, whereas 
y- is the outward normal to the ramp. The SL surface is 
not flat and given by 

y — 4>(x), 4>{x) = S sin kx. (1) 

AL surface is given by 

y = £(x) = <j){x) + h(x), with < h(x) >= h , (2) 

where < h > means average over wavelength in x- 
direction. 

In the experimental data surface velocity of fluid 
is about 3cm/s, and by using Nusselt equation which is 
shown later g , the average water-layer thickness Hq was 
around 0.1 mm. The wavelength A = 27r/fc is about 1 cm 
experimentally. We use the parameter 

fi = kh , (3) 

which is 6 x 10~ 2 for the typical experimental values 
above. In general, having /j, « 1 defines the long- 
wavelength approximation. 

Our key assumptions are as follows. 

• Steady-state(time-independent) flow. 

• /i << 1 (long wavelength approximation) 

• Incompressible fluid. 
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We prefer to use the following dimension-less variables. 
For length, 



y -> y*= y/h . 

The thickness of water layer is thus 

h(x) -> h*(x) = h(x)/h = 1 



h(x) 



(4) 
(5) 



(6) 



And the height of the SL and AL surfaces are 



SL : (f>*(x) = {5 /ho) sinx* = 77 sin x* 



(7) 

AL : £*(x) = h*(x) + <j)*{x) = 1 + h(x) + ?7sinx*(8) 
The characteristic flow velocity in the x-direction is 



Uo 



gh,Q sin 
2v ' 



(9) 



where g is the gravitational acceleration and the viscosity 
v = 1.8 X 10~ 6 m 2 /s at degrees centigrade. This is 
a Nusselt equation, the theoretically predicted velocity 
at AL surface when S — 0, i.e., for a flat SL surface 
with uniform thickness || . This unperturbative solution 
results by equating the gravitational force to the viscous 
force. The velocity distribution is parabolic in y; 



U (2 



y_ 

ho 



(f) 2 )- 
ho 



(10) 



We consider perturbations of this solution. By us- 
ing the above formula for the speed, we relate the 
experimentally-determined flow Q to ho and, surface ve- 
locity Uo- The flow quantity is defined by Q — 2irRUho, 
where U — J Q h ° v x dy = 2E/o/3, the mean speed; R, the 
radius of icicle. In the experiment ||, Matsuda used the 
flow Q = 160ml /hr and width of gutter I — 3cm (/ corre- 
sponds to 2ttR) because this produced the clearest waves. 
This gives Uh = 1.48 x lQ- Q m 2 /s. From U = 2U /3 

and Uo = ^ (by setting 9 = x/2), we get Uo — 
2.4 x lQ- 2 m/s with ho = 0.93 x 10" 4 m. On the other 
hand, his measurement of the surface mean velocity by 
observing the motion of fine grain was Uo — 4 x 10~ 2 m/s 
at 6 = 7r/2. Hence, we assume U = (2.4 ~ 4) x 10 _2 m/s 
with ho = (0.93 ~ 1.21) x 10 _4 to as the experimental 
surface speed and water-layer thickness. 
In the y-direction, characteristic velocity is 



(11) 



We write u as speed in the x-direction, v as that in the 
y-direction, and P as the pressure. Dimensionless speeds 
and pressure are given by 



u/Uo, 
v/Vo, 
P 



pgho sin 9 



(12) 
(13) 

(14) 



Other dimensionless constant parameters are, respec- 
tively, the Reynolds number and the Weber number 



R 
W 



h U 



pyhl' 



(15) 
(16) 



where 7 ~ 7.6 x 10 2 N/m is the surface tension of 
liquid water. Approximate values of these quantities, 
U ~ 3 x 10" 2 to/s, h ~ 10- 4 to, and v ~ 1.8 X 10~ 6 m 2 /s 
predict R ~ 1.5 and W ~ 10 3 . This value of Reynolds 
number indicates laminar flow. The flow components u 
and v satisfy the steady-state Navier Stokes equation for 
incompressible fluids: 



(v . V)v = - 



VP 



(17) 



The incompressibility condition is 



V-v = 0. (18) 
Mass conservation law at the water-air interface requires 



dx u(x, £(#)) 



(19) 



The boundary conditions at the SL surface are zero fluid 
velocity: 



u(x,y = <f>(x))=0, v(x,y = <f>(x)) = 0. 



(20) 



Stress balancing condition on AL surface; that is, free 
surface condition is 



P{in) n i = Pn % 



pv ( dvi + "~ K )n k 



dvk . 
dx 1 



7" 



d 2 H(x) 
dx 2 



(21) 



The index means x 1 = x\ x 2 = y and rii is the normal 
unit vector to the AL surface. P(i n ) is the pressure just 
under the AL surface, and P is the atmosphere pressure. 
For the above equations, we approximated surface ten- 
sion term as 7^" by neglecting the second order term in 
M- 

Now we rewrite the equations in dimensionless form. 
The incompressible fluid condition is 



dx* 



dv* 
9y* 



= 0. 



To hold this condition automatically, we introduce di- 
mensionless stream function by 



dtp 
dy* 



v* = — 



dip 
dx* 



(22) 



In the following, we use the stream function instead of 
the velocity. Also, we drop the * mark on the dimension- 
less quantities. All the quantities are dimensionless till 
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the end of this section. 

Now Navier Stokes equation becomes 

t/jyyyy = fJ,R[tf, <y^ > X yy ~ ^X^yyy] ~ 2fi tj) 

xxyy 

X XX 

] - fj, 1p 

xxxx • 

(23) 

where the indices indicate derivatives with respect to x 
andy. The 4th-order derivatives appearing in l.h.s. come 
from viscous term with taking derivative to cancel out the 
pressure term in equation (fl7|). The pressure is deter- 
mined from Navier Stokes equation by using the stream 
function as follows. 

P * = ~ + Y^vw - §U J y^y ~ Myy) + ( 24 ) 



and the stream function given by 



if, 



1 



(y — rjsinx) 3 + (y — rj since) 2 



(35) 



In our approximation we have only 0-th order terms 
in /i. This is because the first order terms in /x are also 
proportional to the small quantity rj, the amplitude of a 
small ripple, which makes them effectively second order 
quantities. The form of stream function is intuitively 
understood easily. Because it is just the modification of 
unperturbative solution for the velocity to be vanished 
at non flat SL boundary. 



III. THERMAL DIFFUSION PROCESS IN THE 
WATER LAYER 



Py COt ij-Jxyy ^ ^xxx 



t i 2 R 



(25) 



For the stream function we include only 0-th and 1st 
orders in /x, and for pressure, only 0-th order terms are 
kept. 

The Navier-Stokes equation for the stream function 
and the pressure equations become 



fiR[lp v 1pxyy ~ 4>x4>y 



(26) 



11 R 

P x = - + l^i>yyy--^{ll>yi>xy-'4>x'4>yy), (27) 



P, 



cot ( 



(28) 



Next we consider the boundary conditions up to O(fi). 

(29) 



ip x (x,y = (j)) = ip y (x,y = cf)) = 0, 
P(x,y = = t 

lf}yy{x,y = 0=0 



Wo - 

P{z,y = = P ^-z(hxx - V^inx) - fiip xy (30) 

smu 

(31) 
(32) 



h x + r]cosx = -^-(x,y = £), 



where we define Wo = [i 2 W ~ 10° as order 1, because W 
is 7.6 x 10 2 in our case. 

When rj = h = 0, flat laminar flow occurs with the 
solution 

ip = -\y 3 + y 2 , P = P + (l~y)cot9, (33) 

which is easily shown to satisfy Navier Stokes equation 
and all boundary conditions. 

Therefore, we consider the perturbations from this so- 
lution. The precise perturbative calculations are given 
in appendix, and as a result we obtain the height of AL 
surface, 



£(x) = 1 + rjsinx, 



(34) 



A. Basic equation 

Now we consider the thermal diffusion process in the 
fluid that was obtained in the previous section. We make 
the following assumptions. 

1. We use the long wavelength approximation up to 
first order in fj, (/x = kho ~ 6 x 10~ 2 ). 

2. We ignore thermal expansion of the water and thus 
retain incompressible fluid. 

3. Heat transport is through steady-state thermal dif- 
fusion with flow. 

Note that steady-state is valid because the time scale for 
temperature change is much longer than the time scale 
for the ice crystal growth. The heat flow is given by 
Q = —nS/T+(pcT)v, where k is the thermal conductivity 
of water, T the temperature, and c is the specific heat of 
water. The steady-state continuity condition is given by 
dropping the time derivative. 



AT VT = 0, 

D 



(36) 



where D = — is the thermal diffusivity of water and the 

pc J 

incompressibility condition was used. Furthermore, we 
dropped the term for the thermal energy coming from 
energy dissipation of fluid [pi , 



2k dxk dxi 



(37) 



Because this term is much smaller than the other terms. 
In a mean while we use dimensionless parameters (a;*, y*) 
again. From the previous section, 



u = U 



dtp 
dy* ' 



-/j,U 



dtp 



where tp is the dimensionless stream function. Then 
equation ([56]) becomes 

^1 = a \^— - ^L^L] (38) 
dy 2 dy* dx* dx* dy* 
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where we have dropped a /i 2 term. In experimental data, 
D ~ 1.3 X 1CT 7 , /j~6x 10~ 2 , f7 ~ 3 x 1(T 2 , h ~ 1(T 4 . 
This gives 



a = fi- 



h Uo 
D 



1.4. 



In the following sections we drop * mark on the dimen- 
sionless quantities again. 



B. Expansion in powers of y 



We start from equation (|38| ) with stream function fl35[ ) . 
The temperature at SL boundary will be very nearly 
equal to the melting temperature Tm = 273.15.fC at at- 
mospheric pressure. The surface tension for the curva- 
tures in the experiment (Gibbs-Thomson effect ) can alter 
the melting temperature by at most 10~ 6 K, which can 
be neglected. Hence we expand the solution in powers of 
Y = y — r\ sin a;. 



T(x,y) = T M + a 1 {x)Y + a 2 (x)Y 2 + 



(39) 



The left and right sides of (38) become 

oc 

Tyy = T Y y = ^2 k ( k ~ lW(x)Y k - 2 , (40) 



k=2 



dip dT dip dT 
dy dx dx dy 
dip dT dip dT 
dY~dx ' Y ~~dx~ ' Y dY 

{ 2Y-Y 2 )jr^Y k 



k=l 



dx 



(41) 



Setting the left and right sides of (|3q ) equal gives 

oo oo 

k(k-l)a k {x)Y k - 2 = aY(2-Y) 



da k (x) YK 



k=2 k=l 

The coefficients are found recursively. 

a d 



dx 



1n+4 — 



(n + 4)(n + 3) dx 



{2a n+ i — a n }, 



adaAx) 

02=03 = 0, a 4 = — — . 

b dx 



(43) 



All the coefficients arc determined when a\{x) = a(x) is 
known. The first nine coefficients are the followings by 
using the definition D = a-4-. 



a i 

"3 



a, 
0. 
0. 



a 4 

«(-> 

a 7 



Of) 



D 

J"' 
D 

-To 1 
o, 

D 2 

b 2 



126°' 



Because a ~ O(l), we 
second- derivative terms. 
On the SL surface, 



210 

D 2 

a 

1440 



consider only through these 



T(SL) = Tm = const. 

dT 

Q(SL) = \y=r)s\nx- 



— Kd{x). 



On the AL surface, 



T(AL) 
Q(AL) 



T, 



M 



[1 



-k[1 



12 



60 
-D 



-D 



13 



3360 



D 2 ]a(x), 



239 



D 2 ]a{x), 



(44) 
(45) 

(46) 
(47) 



10080 

where Q is the heating resulting from the tempera- 
ture gradients. We have omitted the 0(b 3 ) terms. 
These additional terms are (6.0 x 10~ 5 )D 3 a in ([i"6|), and 
-k(5.2 x 10- 4 )L> 3 a in @. To be comparable with 
0(b 2 ) term, a needs to be about 10 2 . Therefore, this 
approximation is valid when a <C 10 2 . To determine 
a(x), we must consider the temperature and heat flow at 
the AL surface. 



IV. THERMAL DIFFUSION IN AIR 

We consider the thermal diffusion in air to consider 
the temperature and heat flow at the AL surface. We 
note two points here. First, we can not approximate 
the icicle system as the ramp. The picture of ramp is 
a good approximation when we consider inside the thin 
water-layer, but not good for outside. Therefore we treat 
the icicle as cylindrical object, and consider the thermal 
diffusion outside. Second, we can not use the same di- 
mensionless variable as before since our physical space 
is outside. Therefore we use different dimensional coor- 
dinates in this section. The diffusion equation in air is 
given by 

AT = 0. (48) 

Let us write down in cylindrical coordinate. 

d 2 Id Id 2 d 2 
[^ + -^ + ^ + ^]nr,M)=0. (49) 
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We assume axial symmetry, and so we have dT/dO = 0. 
Therefore, we work with 

d 2 Id d 2 
^ + r8-r + ^ T ^ = °- (50) 

Because the surface oscillates in the x-direction, we as- 
sume that the solution has the form, 



T(r, x) — f(r) + g{r) sin(fcx + (/)), 



(51) 



where we have assumed that the icicle is an infinitely long 
column with small surface fluctuations. f(r) satisfies 



, d 2 1 d , „, . 



and g(r) satisfies 



r d 2 Id , , 



(52) 



(53) 



The solution is, 

T(r, x)=A + B log(r/i?) + CK (kr) sm(kx + 0), (54) 

where Kq is the 0-th modified Bessel function, and A, B 
and C are the constants. R is the mean radius of the icicle 
including the thickness of water layer ho- Note that the 
constant C is the order S, because it is introduced from 
the fluctuation on icicle. We define local coordinate y by 



R + y. 



The solution is 



T(x,y) = A + Blog(l + y/R) 

+CK (k(R + y))sm(kx + (j)), (55) 

Note that the ramp system is retained by taking the limit 
R — » oo with fixing B/R and Cexp(— kR) finite. 



Near the AL surface, y << R, we have 

T(x,y) = A+-y+--- 

+ C(K (kR) + K' (kR)ky + • • •) sin(fca; + 4>) 
~ [A + CK (kR) sin(kx + <f>)] 

+ [— + CkK' (kR) sin(kx + 4>)}y 
R 

+ \[~ + Ck 2 K' Q \kR) sm(kx + cf>)]y 2 + ■ ■ ■ ,(56) 

where K 1 and K" indicate derivative of K with respect 
to its argument. We define the mean growth rate of icicle 
V^by 

\r - K o ^ dT n B 

V --T < ^ > y= 0= -TR^ (57) 

where kq is the thermal conductivity of air, and < > 
means spatial average over x. So we obtain 

LRV 



B 



(58) 



At y = Ss'mkx (AL surface), the temperature and heat 
flow are, 

T(AL) = [A + CK (kR)sin{kx + (/))} 
LV 

+ [ h CkK' Q (kR) sin(kx + 4>)]5 sin kx 

Kq 

LV 



= A + [CK (kR) cos 4> 5} sin kx 

+ CK~o(kR) sin (j) cos kx, 



(59) 



Q(AL) = LV — [KoCkK' (kR) cos < 



LVS 
R J 



sin kx 



KoCkK' n (kR) sin 4> cos kx. 



(60) 



Because the constant C is the order <5, we have dropped 
5C terms, and we keep terms up to first order in 5. 



GROWTH RATE 



T(x, y)= A + B'y + C cxp(-ky) sm(kx + <p), 

where B' and C are other constants. Hereafter we 
work with finite R case, because the ramp case is always 
retained by taking the limit as above. The reader might 
think the appearance of logarithmic term strange, since 
it diverges at large r. But the appearance of such a term 
is natural for infinitely long axially symmetric source. As 
real icicles have finite length, this solution is valid only 
close to the icicle; far from the icicle, the icicle acts like 
a point source of heat and we must mach the near-icicle 
and far-icicle solutions at r ~ L where L is the length 
of the icicle. This matching includes the temperature at 
infinity and partly determines the coefficients A, B and C; 
in addition, the mean growth rate of the icicle radius also 
determine coefficient B, which includes the information 
of the temperature at infinity. 



Two boundary conditions apply to the AL surface: 
continuous temperature and continuous heat flow across 
the water-air boundary. Comparing these two sets of 
equations ([hj]) and (^9|) , (^) and (p0|), the solution re- 
quires that a(x) = E + F sin kx + G cos kx, where E, F, 
and G are constants with dimension of temperature. 

Then, in dimensional units, ( |46| ) and ( |47| ) are 



T(AL) = 



M 



E 



[G- 



13a 2 
3360 
13a 2 



F] sin kx 



—F 

60 3360 



G]coskx, (61) 



Q(AL) 



-E 



ho 



[G- 



-[F- 
i 

5a 
h 12- 



5a „ 239a 2 



-G 



F 



12 10080 
239a 2 



F] sin kx 



10080 



G]cosfca;. (62) 
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By comparing to (^9|) and (|60|), we have six equations 
with six unknowns, ( i.e., A,C,(/>,E,F, and G). A and E 
can be written immediately: 



A = T, 



LVhn 



M 



E = — 



LVhr, 



(63) 



(64) 



The other four quantities are determined by the following 
equations. 



1 



_ 13a 7a 
3360 60 „ 

7a i _ 13a 2 

60 1 3360 




CK (kR) COS(f>- ^ 

CK (kR) sin^ ° 



(65) 



1 - 



239ct 
10080 
5a 
12 



1 



_ 5a 
239a 2 




10080 

( tfixK' {kR)C cos 4>+ L ^ 1 
ffiK^(kR)Csia(f> 



Solving above equations, we obtain 
LV5K' Q 1 



/'- 



239 2 
10080" 



G = -n 



kK 1 
LV8K' n 



1272 2 
10080" 



( 239 y 4 
V 10080 / " 



kK 1 



2272_q,2 , fJ39_)2 Q 4 : 
10080" ' V 10080/ " 



(66) 

(67) 
(68) 



where we have neglected second order terms in /i and we 
neglected a term proportional to ho/R(~ 10~ 3 ). 

At the SL surface, the growth rate of ice is given by 

= Q^SX) = v _ k , Fsinkx + Gcoskx \ ( 69 ) 
L Lho 

The form of growth rate is different to usual MS the- 
ory. In MS theory, v(x) = V + /sin foe for the surface 
y = 5 sin kx. But now we have another term cos kx. To 
understand its physical meaning, we write the relative 
growth rate v s as the growth rate in a reference frame 
moving with velocity V (v s =v—V). 



v s (x) = f sin kx — g cos kx, 

LriQ Lho 
The steady-state condition means, 

dy s (x,t) 



(70) 
(71) 



v s (x) 



dt 



where y s is the hight of SL surface in reference frame. 
From the steady-state condition, the time scale for 




hf.Ud 



FIG. 4: The dimensionless amplification factor 
the dimensionless wave number a. Laplace instability in air 
and a GT like effect due to hydrodynamics amplifies pertur- 
bations that have a wavelength given by an a close to 2.2. 



growth of fluctuation is very long compared to our ob- 
serving time scale. By solving the equation: 



v s (x) 
with 

we obtain 



dy s {x,t) 
dt 

y s (x,t = 0) 



t=o= / sin kx — g cos kx, (72) 
Ssinkx, (73) 



y s (x, t) = S(t) sin(A:a; — ojt), 



with relations 



5 = f, u> = g/S. 



(74) 



(75) 



The essential point is that the fluctuation is not only 
growing up, but also traveling downwards (/ > 0, g > 0). 

Therefore the amplification factor is determined from 
F. By using the relation — K = K\ > 0, we have 



LhnS 



:F = Vk- 



gl(fcg) /1 
K (kR) y 1 



239 
10080 



a 2 ) 



1 



1272 



+ (^39_)2 a 4 
' V 10080^ " 



(76) 



Note that #4^t 

K a (kR) 



1 + 2^7j- In the case of kR, » 1/2, 
meaning a thick icicle, we obtain 



1 



239 2 
10080" 



1272 



, / 239 p 4 
~ V 10080 I " 



(77) 



Because a oc k, this form is similar to the amplifica- 
tion factor given by Mullins-Sekerka theory. The ampli- 
fication factor increases proportional to wave number by 
thermal diffusion in air as expected for a Laplace insta- 
bility, and it decays by interaction with fluid (aterms). 
The thermal diffusion in thin water flow works just like 
Gibbs-Thomson effect, since the fluid changes tempera- 
ture distribution to be uniform and inhibit Laplace insta- 
bility. From these two effects, we have maximum value 
for 5/6. 

The maximum amplification factor occurs when a = 
2.2, which determines a preferred wavelength. 



2.2. 
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FIG. 5: The dimensionless speed that the ripples travel 
down an icicle versus a ~ Uq/X, where V is the growth rate 
of icicle. 



per hour in experiment . Because water flow carries heat 
flow downwards, this kind of motion seems to be natural. 
Actually the speed w = when Uo = (a = 0). Then 
we can imagine that such speed will be increasing func- 
tion of fluid velocity Uo- But this is not true. For (solid) 
surface fluctuation to change its form, it is necessary to 
emit latent heat non uniformly. Like the amplification 
(growing) factor becomes maximum at the specific value 
of a, the traveling speed may have such dependence of a. 
Qualitatively this is due to the fact that fluid not only to 
carry heat downwards (small a region in fig. 5), but also 
make the temperature field uniform and thus suppresses 
the heat diffusion (large a region in fig. 5). 



By using D = 1.3 x 10 7 m 2 /s with experimental data: 
Uo = (2.4 - 4) x lQ- 2 m /s, h Q = (0.93 - 1.21) x 10" 4 m, 
we have 



2n h 2 U 

- = Z7T 



5mm ~ 13mm, (78) 



which well agrees with experimental value of 8 mm. ||] 
Also, in agreement with observations, this value does not 
directly depend on external temperature. 

Using the constraints between Uo, ho, and flow quan- 
tity Q as, 



Uo = 9 -^, Q^RUoho. 
2v 3 

Then Q is related to A max as 



(79) 



A„ 



1 / v M/3/ 3( 3 



(_)l/3(^.)4/3 



Dctmax gn 2i? 



where = ir/2 is assumed. 

Therefore if Q is proportional to R for usual icicles, its 
wave length is uniquely determined. But for thin icicles 
with R < A/(47r), we should include the term that 
we have neglected. 

Next we consider the travel of fluctuations along the 
icicle. The traveling phase velocity w is the following 
function of a. 

-! = lf^ = §f *«>~V««), (80) 
where /3(a) is defined by 

5 



/3(a) 



12 



a 



1 



2272. 2,/^39_p 4' 
10080 ~ \ 10080' 



(81) 



and its form is given in Figure 5. 

At the 8-mm wavelength, fluctuation travels down- 
wards very slowly with speed w ~ 0.5U, where V is the 
mean growth rate (V = R) . For the ambient air temper- 
ature of -8 degrees centigrade, this speed is about 1 mm 



VI. CONCLUSIONS 

We have shown that an icicle covered with thin water 
flow makes wave-like undulations on the ice during so- 
lidification, and the preferred wavelength is determined 
by a MS-like theory. Thermal diffusion in air makes 
wavelength shorter, i.e. the amplification factor becomes 
larger for shorter wavelengths. On the other hand, ther- 
mal diffusion in thin water flow makes the wavelength 
larger, i.e. the amplification factor becomes smaller for 
shorter wavelengths. From these two effects, a specific 
wavelength emerges with a maximum amplification fac- 
tor of fluctuation. The thermal diffusion in the thin wa- 
ter flow works just like the Gibbs-Thomson effect because 
the water flow makes the temperature distribution more 
uniform and thus inhibits the Laplace instability. This 
then is one of our main results. 

Our A max depends on 9 and flow quantity Q = 
2-xRhoU — lh U. (8 dependence will be observed for 
flow on inclined gutter experiment: Fig2. For an icicle, 
6 = tt/2 should be used.) 



1 



Da r , 



gn sin < 



]) ( 2R ) , 



(82) 



where we have used the relative Nusselt equation M, 



ho = [ 



vQ 



11/3 



guR sin 9 



(83) 



The 9 dependence of the wavelength for the ramp case 
is discussed by Matsuda 0] experimentally as 



A 



8.2 



sin ' 



(84) 



with 7 = 0.6 ~ 1. However, the number of plotting 
data is few, and thus we can not give a definite result 
experimentally. 7 is 0.3 in our theory. The difference 
exists, but not a big disagreement qualitatively. 

In nature, all icicles have their own flow rates Q, but 
almost all have the ripples with the same wavelength. 
This is explained by our analysis because the wavelength 
depends on ratio of Q to R, but not only on Q itself. 
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It is then natural to assume that in nature, Q is propor- 
tional to R . The selected wavelength of our analysis does 
not depend explicitly on external temperature. Although 
the mean growth rate V and also the amplification fac- 
tor increase with a decrease of temperature, the selected 
wavelength with the maximum amplification factor is in- 
dependent of icicle growth rate (j77|). 

We have also shown that surface ripples are expected 
to travel downwards during icicle growth with a speed 
of 0.5 times the average normal growth rate of the ice. 
This should be checked by experiment. Our theory may 
be used to map the waves around mineral stalagmite by 
changing the diffusion equation from temperature field to 
a solute density field. Furthermore, our diffusion equa- 
tion in fluid is mathematically similar to the Schrodingcr 
equation for a harmonic oscillator having the complex 
valued potential. Therefore, the algebraic method may 
be possible to be used for analysis. This point will be 
discussed elsewhere. 
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VII. APPENDIX 

We define the fluctuation fields tp and P by the follow- 
ing relation: 



^ = -^y 3 + y 2 + ${x,y), 

P = (1 - y) cot 9 + P + P. 



(85) 
(86) 



Then for the fluctuation fields, we have following equa- 
tions: 



R w 2 

-2-K-f 

-{-2y + 2 + tjj yy )iP x }, 



1 R 

Px = 1^4>yyy - l^[(-y 2 + 2y + *i>y)i>xy 



P y = 0. 



The boundary conditions are 
MSL) = 0, 



tpy (SL) = rj 2 sin 2 x — 2rj sin x, 
P(AL) 



(87) 

(88) 
(89) 



(90) 
(91) 



\hxx - V sin x) 



+ (h + r) sin x) cot 6, (92) 
^yy(AL) = 2(h + rj sin x), (93) 
jj) x (AL) = (h x + rjcosx) x [(1 + r/sinx + h) 2 

-2(1 + rj sin x + h) -4> y (AL)], (94) 

where the AL surface is determined by y — 1+rjsm.x + h, 
and the SL surface is determined by y = rj sin x. 

To solve the above equations, we assume the solutions 
have the form, 



P 

ryyy 



p(0) + . . . i 

0, 



(95) 
(96) 
(97) 



and neglect higher order in fi. The solution for the 
0-th ordered stream function and pressure are 

= (h + rj sin x)y 2 - [rj 2 sin 2 x + 2(1 + h)r] sin x]y 
1 



+ [{l + h)rj 2 sin 2 x + -rj 3 sh 



(98) 



p(0) 



Wq ~ 

: [hxx — V sin x) + (h + rj sin x) cot 9. (99) 



sin f 



To obtain the lst-order stream function, we put above 
solution into the following equation. 



(i) - 



2Wn /7 



an t 



(hxxx — V cos x) + 2(h x + rj cos x) cot I 



+R[(-y 2 + 2y + $°>)$°> - (-2y + 2 + ^^frOO) 



which is obtained by ( |8q ) with the help of (|99[). This 
expression is consistent with (jS7|). After some integra- 
tions with the boundary conditions, we obtain 

= — {I + h)h x y° --rj{l + h)h x sin xy* 

R ~ ~ 

+ -rf (1 + h)h x sin 2 x y 3 

1 W ,~ 

+ ol — • 7i i^xxx V cos x) + (h x + rj cos x) 
6 Sin a 

x cot0] y 3 + \c x (x)y 2 + C 2 (x)y + C 3 (x)(101) 

C\,C2,Cs are determined by the boundary conditions. 
The result are 



Ci = -2(1 + h + rj sin x)[R(l + h)h x 

x {-(1 + h + r/sinx) 2 — (1 + h)rjsinx} 

(hxxx - ?7C0SX) 



sin 9 

+ (h x + rj cos x) cot 9], 
C* 2 = —R(l + h)h x rj 4 sin 4 x 



(102) 



sin 9 



W - 

+ \-^-7,(h X xx - ijcosx) - (h x + Tj cos x) 
sm u 
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x cot 9]rj 2 sin 2 x — Ci(x)r] sinx, 
dC 3 R r 
dx 

R 



(103) 



'■[{I + h)h x ] x rf sin 5 a; 



~(1 + h)h x rf sin 4 x cosx 

1* r W ° (I ^ ■ \ 
3 smO 

— — jysinx) cot 0] sin 3 a; 

77 2 

— — d x Cx{x) ■ sin 2 x — r\d x C2{x) • sinx. (104) 
Now we put ip = ip^ + fi'ip^ into (pi). We now have, 



2(l + h) 2 h x + ^ x 1 \AL) 
+ (h x + rj cos x)^{AL)] = 0. 



(105) 



We expand the fluctuation of water-layer thickness as 
h = U a) (106) 

and by putting it into above equation, we have h x = 0, 
and so we can use 



hf® = 0. 



(107) 



This can be done by redefining of h$ after subtracting 
a constant. Then m 1 ) is determined by 

hi 1] = -\$£HAL) + r ? cos^W(A£)] s=0 , (108) 

From the fact that h starts from 0(p,) in its expansion, 
we can determine ip^ as , 



(i) 



W 1 a 

(-r-5 + cot 0) [—77 cos x ■ y 
smv 3 



•q cos x(l + r/ sin x)y 2 + {rj 3 sin 2 x cos 2; 

+ 2?7 2 sinx ■ cos x}y] + C 3 , (109) 



9xC 3 = 77 sina;(— 



VK 



sin 6* 



cot ( 



x [—r; sin x — 77 sinx + 3 sin x — 2j. (HO) 



From the above expression and (108), we obtain the 
following simple relation: 



h x 1} = ~(^ + cot 8) sinx. 
3 sinfl 



Furthermore, 



= !^(J^_ + cot e \ cos ^ 

3 sinw 



The height of the AL surface is, 



(111) 



(112) 



£(x) = 1 + rysinx, (H3) 



whereas the stream function is 



V> = -\i/ + y 2 + ^\h) + ^ {1) , 

= (y — rfsinx) 3 + (y — r/sinz) 2 . (114) 
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The amplification rate of a wave (fluctuation) depends on 
its wavelength. Only the wave with the largest amplifica- 
tion rate can survive during the development of the fluc- 
tuations, other modes will be observed as lower-amplitude 
noise. The constant rate of 5/5 means that the growth is 



exponential in time; hence, small differences of this con- 
stant results in large differences in fluctuation amplitude 
at long times. In this sense, a mode distribution like that 
in Fig. 4 will be hidden in the external noise and thus not 
observed in experiments. This situation is the same as in 
the original MS theory. 



